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Abstract. A review of current state of understanding of dielectric mixture 
properties and approaches to use numerical calculations for their modeling are 
presented. It is shown that interfacial polarization can yield different non- 
Debye dielectric responses depending on the properties of the constituents, their 
concentrations and geometrical arrangements. Future challenges on the subject 
are also discussed. 



1. Introduction 

Many materials used in our daily life are composites, which are often made up 
of at least two constituents or phases. The outstanding mechanical properties of 
many composites, and especially the unique combination of low density with high 
strength and stiffness have led not only to extensive research but also to a highly 
developed technologies |[ ||]. In comparison, relatively little attention was given 
to their other physical properties, which in parallel affected their use in electrical 
appHcations |, |, |, 0]. The main advantage of composites is the ability to tailor 
materials for special purposes. Their applications are evolving day by day through 
the developments which lead to better precision in their manufacturing. 

Designing of composite materials for electrical applications with classical trial 
and error approach requires a lot of time and money. However, by using computers, 
modeling desired properties of insulation systems can be estimated and corrected. The 
electrical properties of the system, i.e., its conductivity and dielectric permittivity, 
are influenced by the properties of the constituents, interaction between them 
and geometrical configuration. One should not forget that insulating materials or 
dielectrics, as Faraday called them show various properties at different voltages, 
temperatures, frequencies, moisture content and mechanical stresses. These should be 
considered in the design as well as in the diagnostics. 

From the early days of the electromagnetic field theory, predicting and calculating 
the dielectric properties of mixtures has been a challenging problem of both theoretical 
and practical importance @, ||, ||, |6[ 11, ^ In composites the polarization 



of charges due to the differences between the electrical properties of the facing 
constituents plays an important role. Maxwell |jl^ was the first to notice this 
phenomenon when he considered a binary layered structure and expressed the effective 
dielectric properties of the composite. When structures other than the layered 
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arrangement of phases were studied |T^, ||, ^, |T^, it was realised that the 
problem was not trivial, and even the geometrical shapes 18, |l^, |l] and 
arrangements of inclusions ^ ^ ^ played an important role 

in the electrical properties of composites. 

The aim of this article is to give a brief review on dielectric mixture systems 
studied in the literature and to attract attention to the behavior of interfacial 
polarization in frequency dependent fields. First, dielectrics and polarization 
phenomenon are presented in the next section. Later sections summarize different 
theoretical studies available in the literature as well as numerical simulations 
performed by the authors on ideal structures. The simulations presented were in 
two dimensions and were designed such that the chosen material parameters have 
yielded dielectric relaxations at frequencies lower than 1 MHz. Simulations of three 
dimensional structures are so far very rear (29[ |3^, |3^, ^ and there is still a need 
for developments. 

2. Dielectrics 

2.1. Polarization in dielectrics 

The interaction between electromagnetics fields and matter is described by Maxwell 
equations. Vectors expressing the electrical components, dielectric displacement D 
and electric field E of the electromagnetic phenomena in dielectrics are interrelated. 

D=eoE + P (1) 

The constant Sq is the permittivity of the free space, l/367r nF/m. The polarization 
vector, P, can also be written as follows, 

V-P = -p (2) 

where p is the charge density. The relation between the dielectric displacement, D, and 
the applied field, E, is often linear and can be expressed with a simple proportionality 
constant, e, 

D =eeoE (3) 

The constant e is called the relative permittivity, which describes the dielectric 
properties of medium. When the polarization, P, is taken into account (inserting 
Eq. in Eq. (^), it is also proportional to the field, E 

P=X£oE = £o(£-l)E (4) 

The quantity x is called the polarization coefRcient of the substance, or its dielectric 
susceptibility. Observe that in free space, s = 1, there is no polarization and x = 0- 

The polarization in materials can be due to several mechanisms; electronic, 
ionic (molecular), atomic, dipolar (orientational), and interfacial polarizations [ p^ , 

|36| , |37[ |. Furthermore, hopping of charge carries between localized sites also 
creates polarization |3^, ^ The first three polarization mechanisms are 

much quicker than the others. For this reason, depending on the time scale, the 
fast polarizations can together be considered as an instantaneous polarization, Pq. 
Moreover, the polarization is to be finite at longer times, P{t oo) = Pg, since 
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Figure 1. Time dependence of polarization, P, when a constant electric field is 
applied at t = t;. 



infinite amount of electrical energy can not be stored in a dielectric (or a capacitor). 
Fig. |l| shows the polarization for a medium when a constant-step electric field, E = i?on 
with n being the unit vector, is applied at t = ^. The permittivity values at t oo 
and t — ^ are defined as Eg and £oo which are the static and instantaneous (high 
frequency) dielectric constants, respectively. The time dependence of polarization is 
expressed as following when the electric history of the material is known 

P(i) = £o(£oo - 1) E(0 + eo f Ht - m)d^ (5) 

where f(i) is called the dielectric response function and it is an intrinsic material 
property. For harmonic fields Eq. (|^) and Eq. (|l|) yield to the complex dielectric 
permittivity, e, as follows. 

/•oo 

e{uj) =£00 + / f (i) cxp{—iLut)dt (6) 
Ja 

Therefore, the Fourier transform of the response function f (i) is defined as the complex 
dielectric susceptibility, 

X(w) = e-£oo 

f°° (7) 
X{uj) = x' - «X" - / f (0 exp(-^c^^)d^ 

The real and imaginary parts of the dielectric susceptibility are coupled to each other 
by Kramer-Kronig relations |^ . 

2.2. Interfacial polarization 

When two media are put into contact (forming an interface) and an electric field is 
applied, charge polarization occurs at the interface due to the differences between the 
ratios of the electrical properties (conductivity and permittivity) . This phenomena was 
first studied by Maxwell 1 4^] who considered two phase system with one of the phases 
insulating, ai = 0. The system was a layered strucure and the effective dielectric 
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permittivity of the mixture was expressed as a function of intrinsic electrical properties 
of the phases, (Ti_2 and ei^2 = £001 2 ^^'^ their volume fraction q. The complex dielectric 
permittivity of the phases, e, in that case is 

ei{uj) = £1 

^ ^ ^2 (8) 

and the effective complex permittivity of the mixture can be expressed as 

~ / ^ 

ee(^) = ee + — (9) 

1 + tUJT 

where Ae is the dielectric strength of the interfacial polarization and t is the 
relaxation time of the polarization. This expression is in the form of Debye relaxation 
function [||. 

{l + q)e2 (l + g)cT2 ^ ' 

and ee is the effective dielectric constant at optical (high) frequencies 

£162(1 + 2g) 



q{ei + £2) + £2 



(11) 



A more general case is that both media have complex dielectric permittivities, Si and 
£2, in this case the effective dielectric permittivity is expressed as |l^, 

£e(^) = ~ ,~ (12) 

qei + (1 - q)S2 

The polarized charge density in the interface p can be calculated from the Maxwell 
equations, and is expressed as 



p{uj) oc 



ei£2 - ^2£l 



qei + (1 - q)£2 



(13) 



The approach of Maxwell (two layer structure) is trivial and yields the only fully valid 
analytical formulas available in the dielectric mixture theory for any composition of 
constituents. It can easily be extented to multi-layer systems. 

2.3. Electrical conduction 

In reality, there are no perfect dielectrics, where conductivity is not present. The 
resulting total current density flowing through a dielectric when a step electric field is 
applied, can be written as, 

=eo[eoo(5(0+f(0]E(0+^E(i) 

The two asymptotic parts of the current density, j, are the instantaneous current 
density due to the capacitive component, £o£ao6{t), and the dc conduction current 
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Figure 2. Charge transport between hopping sites. 



density due to the conductivity, cr, of the material, respectively. The current density 
due to the materials polarization is given by the response function, f(i). 

The conduction in materials may be real dc conduction which can be divided into 
two classes: band conduction ^ and dc hopping conduction The band 

conduction is present in the absence of defects and impurities, and it is led by the band 
structure of the material. The latter, dc hopping conduction, takes place via defects 
or impurities which form potential wells (traps or localized states) that are favorable 
for charge carriers (electrons, holes and ions) to hop. It is therefore a phenomenon in 
disordered materials. 

The first approach to hopping conductivity was described by Mott and Davis . 
They had improved the work of Anderson ji^, ^ on localization (trapping) of 
electrons. The hopping conduction model of Mott and Davis showed that, at a finite 
temperature, electrons are able to tunnel between localized states as in Fig. |2| (for 
further reference see [^6| ). The energy and position distributions of the hopping sites 
decide the transition rate, pi-2 @, 



pi_2 ~ exp 



A 

-2aR 



(15) 



where A is a function of the energies of the states -Ei,2, and /cbT is the thermal energy; 
/cb is the Boltzmann constant, /cb = 86.174 |j.eV/K. Note that the charge carriers 
can in principle tunnel to any state depending on the product of inverse localization 
length a and hopping distance R. Although, the distance R in Fig. |^ is between two 
neighboring sites, the model uses the mean hopping distance. For a disordered system 
in zero electric field condition the rate of change of charges moving from site 1 to site 
2 is equal to the rate of change in charges moving from site 2 to site 1. When an 
electric field is applied, the rate of change will be biased to one side depending on the 
polarity of the applied field. 

The localized states are disordered in insulators and amorphous semiconductors, 
when their spatial and energy distributions are taken into consideration. Therefore, 
perturbations in energy of the localized states activate tunneling between the sites. 
Perturbations can be in the form of temperature, pressure or applied electric field. 
If the applied voltage is assumed to be alternating, the frequency of the field affects 
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the hopping rate as weU. In addition, at a particular frequency, the charge carriers 
can hop between two locahzed sites back and forth which is observed as if dipoles 
were present p^ , [sH - Therefore, at low frequencies or higher temperatures, it is not 
easy to separate conduction and polarization mechanisms in insulators. The frequency 
dependent complex conductivity a is given as |l| , 

a{uj) = teoeiuj) 

where is the hopping conductivity. The low frequency dispersion [p4[ or hopping 
conduction contribution, second term on the right-hand-side of Eq. can also be 
included in the complex dielectric susceptibility x 

Xiu:) = —1- (17) 

where and 7 are constant parameters, < 7 < 1- The relationship between n and 7 
in Eq. ^ and @ is 7 = 1 



71. 



2.4- Dielectric loss and relaxation in materials 

2.4-1- Single relaxation time Debye considered noninteracting dipoles in a 
viscous medium. He assumed a response function with a single relaxation time, td, 

fW = — cxpf-— ) (18) 



The time constant td in the equation is called the Debye dielectric relaxation time, 
or simply the relaxation time. The complex dielectric susceptibility, Xj &s a function 
of angular frequency, w, after using the Fourier transform of the response function in 
Eq. 1^ is 

Ae 

X('^) = 



Ae AetuTD ^ ^ 



'l-hw^r^ ^l + w^T^ 



'D ^-TLU Ijj 



However, at low frequencies, losses due to the dc conduction are present, then, the 
susceptibility contains additional term. 



x(^) 



Ae 



" 1 + lu-^tI 



AeuJTD cr 



(20) 



It is possible to distinguish the two independent processes, dipolar relaxation and 
dc conduction using the Kramer-Kronig relation p| |52| , |5^ or using curve fitting 
algorithms based on complex non-linear least squares 



2-4-2- Distribution of relaxation times From the experimental point of view, the 
dielectric response of solid materials, except ferroelectrics, does not often show Debye 
characteristics, and their relaxation characteristics are stretched out in the time 
domain, and similarly in the frequency domain. For this reason, they are rarely 
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modeled by distribution of relaxation times p^ , [59| , pO| |6l| , |6^ or often by a non-Debye 
empirical response ||, ||, ||, |6|, |6^, |67[ g |69[ 0^ ^ ^, ^, ^, 0, ^ . 

The former model considers a distribution function, J-(t), for the relaxation time, r. 
Consequently, Eqs. O and p^, then, become. 



f (0 ^ •^(^) exp ^ j dr 



(21) 



1 + {i^t) 



dr 



^(r) can be regarded as the fraction of polarization processes with relaxation times 
between log r and log r + d log r and it satisfies the condition 



JF(logr)dlogr — 1 



(22) 



There have been several methods proposed to obtain J-{t) ^ ^ ^ 
|6l| , |6^ , ^ |77t , and, the problem is considered ill-mannered. However, using a least- 
squares-algorithm with constraints 1 82 and applying the Monte Carlo technique | ^ , 

a new procedure has been presented |5|]. In the procedure, a least-squares 
method with box-constraints was used and with the application of the Monte 
Carlo algorithm a continuous time-axis was constructed. The survived time constants 
from the least squares yielded the distribution of relaxation times. 

Since it is hard to obtain the actual distribution of relaxation times, empirical 
formulas in the frequency domain based on non-Debye relaxations have been presented 
in the literature |7^, |8^, |8£ , 87, 77 . One of the most frequenctly used can be preseneted 
as: 



[1 + {tU}T)"]f^ 



(23) 



This equation is called the Havriliak-Negami formula |8^ , and a and /3 are empirical 
parameters depending on the shape of the response, (0 < a < 1 A a/3 < 1) | ]67[ |. 
When, a = /3 = 1, we have the classical Debye process. The other famous 
emprirical response functions with distributed relaxation times can also be obtained 
from Eq. (^3|), i.e., for a — 1 and (3^1 the Cole-Davidson and for a ^ 1 and 
P = 1 the Cole-Cole |7^ empirical formulas. In Fig. || different dielectric functions 
are presented as Cole-Cole plots {x"{^) versus x'{^))- 

When the conduction processes are taken into account, the complex dielectric 
susceptibility of a material, x, can therefore be expressed in an empirical form as 



(24) 



where summation over j represents different dielectric relaxation processes within the 
material, or in a more general functional form as 



x(w) 



1 + (iuit) 



dr 



(25) 
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Figure 3. Normalized dielectric relaxation of different functions, Debye single 

relaxation ( ), Cole-Cole with a = 0.6 and a = 0.8 ( ), Davidson-Cole 

with /3 = 0.6 and /3 = 0.8 ( ) and Havriliak-Negami with {a, 0) = (0.6, 0.8) 

and (q,/3) = (0.8,0.6) ( ). 



3. Studies of dielectric mixtures 



3.1. Experimental studies oj dielectric systems 

Techniques widely used in different branches of science for studying dielectric 
properties of materials are the dielectric spectroscopies in time and frequency domains. 
They provide informations on structural and dynamical properties of considered 
systems. There exits in the literature vast of articles on the applications of the 
dielectric spectroscopy to mixtures. A broad range of frequencies, can be accessed. 
The information obtained about the relaxation rates and time constants t (or the 
dynamics) of a system ranges approximately between |j.Hz to PHz in the frequency 
domain and between Ms to fs in the time domain, respectively. The polarization 
mechanisms, such as, electronic, atomic, molecular etc., in the materials can be 
therefore studied. Moreover, the dielectric spectroscopy is an effective tool to 
investigate not only pure form of material states on any scale (micro-, meso- and 
macroscopic), but their interactions as well. 

For insulation applications polymer matrix has been filled with inorganic powders, 
for example silica and aluminum hydroxide, to enhance the mechanical and electrical 
properties. Similar systems can be found elsewhere |8p, |89|, |9C|, |9ll P2l PSi P4|, PSi Pq, 



H, m H, ^ |lO^, [0|, |To|, |To|. Previously it was noted that the dielectric 
permittity of a system is dependent on moisture. In some cases, such as insulation 
of power equipment and components, the moisture content in the system affects not 
only the overall insulation property but the lifetime of the system as well. Therefore, 
measurements of dielectric response can provide information on the moisture content of 
insulating materials |lO|, [O^ |lO|, |lO|, |llO[ |m[ |lT|, [l4[ |lT|| . There are special 
cases where composites have been used for current limiting applications [116, 117, 
lT8| , plgj , [T20I p^ . In those cases knowledge of the dielectric properties is essential. 



In other research fields than electrical insulation, the water uptake o f material s and 
thei r die l ectri c properties have been also studied |l2 
1291 , |l30|, |l3l|. Liquid-sohd mixtures jl3^ 



123, 124, 125, 126, 127, 128 



media il33| , |l3^ , [l35l |7[ H H^, HI, |l2^ 



such as, brine-saturated rocks and porous 



or liquid crystal in polymer matrices |141, 142, 143, 144, 145 1 have been examined 
using dielectric spectroscopy technique. One can also include the mixture examples 



137 , s uspensions in so lutions 1 138 , 139 , 140 1 
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in biological systems where the dielectric spectroscopy has been one of the tools to 
obtain information about the structural and dynamical behavior [146, Jl47, fol |l4^ 



Similarly, work on liquid-liquid mixtures can be found in the literature |149 



153, 154, 155, 156, 157, 158, 159, UM, 161, 162, 163, 164, 165|. The dielectric 



15C, 151 



spectroscopy technique has been found powerful for investigating polymer-polymer 
mixtures |l6|, |l6|, [l69[ |l7^, ^7^ |l7|, |l7|, p74| p75|, p76[ Recently, dielectric 
behavior of confined polymers has got attention [177, |178|7l79, |90|, |18C([. In these 



cases, guest molecules are trapped in regular matrices of polymers or on substrates 
Such structures can also be regarded as composites. 

Charge relaxa tion phen omena in mixtures compose d of ceramics [ 181 
metal-oxides [185, |l86] , ioni c-con duct ors [ 187 , [l88| , 

18£ , 182 , 182 , 184 [ . In the case of conductive coating and 



184 



studied [|185|, |181 



182, 183 



189 , 19C[ have been also 



elec t rom a gnet ic shielding applications carbon-b lack [ 191 , 192 , 168| , 193 , 194 , 195 , 196 



197, 198, 199 and metallic particles [EOOL 201, M have been used as heterogeneities 



in the matrix medium. For photonic applications, metal-dielectric mixtures have been 
widely applied (20| |20|, 1^, |o|, |06[ |207[ |20|, |o|, |l^, |ll[ |212| to obtain photon 
localization. 

Finally, dielectric spectroscopy or similar techniques has been used to sense the 
composition of structures 



213i pMl [2T5I |21^ , |T^, m, Im], m or even buried 



objects [220, 221, 222, 223, 224[. Experimental studies have shown that mixtures 



have enormous potential for special applications. Therefore, numerical simulations of 
mixtures can lead not only better understanding of the physics of dielectrics but also 
improvement of designing of tailored materials without going to expensive attempts. 

3.2. Different analytical approaches 

Predicting the dielectric properties of mixtures has been a challenging problem when 
structures other than the layered arrangements were considered. These different 
approaches are reviewed below. 



3.2.1. The simplest approach Gladstone and Dale [225| expressed a formula for the 
effective electrical properties of mixtures which was proportional and linear to the 
concentrations and dielectric permittivities of constituents 



Ee = (1 - q)ei + <Z£2 (26) 

In reallity, all the expressions that are presented later should yield the same value as 
Eq. ( p6| ) for very low additive concentration {q 
premittivities is not too large. 



0) [ 226 [ , if the ratio of the complex 



3.2.2. Mean field theories The mean field theory approaches suppose that one of the 
phases, an inclusion phase, is inside a matrix phase, and both phases are embedded 
in an effective medium. Wagner [p^ was the first to use this approach, in which in 
the same way as Maxwell |lj], he assumed a structure composed of an insulating 
host medium, ei = ei, with dispersed conducting spherical particles with complex 
permittivity, S2 — ^2 + ij2/{i£olo). The volume fraction of the spherical inclusions 
was q. Later, Sillars [|l^ and Fricke [227[ expanded the con siderations for conducting 
ellipsoidal particles dispersed in an insulating medium [228|. In this case shape of the 
inclusions was accounted for by the parameter n in the complex dielectric permittivity 
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of the mixture. The shape factor n is a function of the eUipsoidal axes and of the 
desired direction lITsIl: 



n, = 



XiX2X3£i 



(27) 



where the subscript i is the desired axis-direction, i — {1, 2, 3}. The vectors Xi^2,3 are 
the radii of the eUipsoid in the arbitrary directions 1, 2 and 3, and Ci is 



dc 



(^? + C)^(x? + C)(a;i + C)(.^l + C) 



1/2 



(28) 



where C is the integration variable. The value of Ci is between 1 and infinity, 
1 < Ci < oo. For needle-like shapes parallel to the field direction, n-values are 
larger than 3. For spherical inclusions, the shape factor, rii, is 3, and for cylindrical 
inclusions if the cylinder axis is perpendicular to the applied field, the shape factor is 
2, and if large disks perpendicular to the applied field are considered, the shape factor 
is 1. As a matter of fact, this ellipsoidal model for the dispersed particles represents 
an infinite variety of shapes which describes real composite materials well. Steeman 
and his collaborators ^ have improved the model of interfacial polarization of 
two phase systems by considering the case when the host medium was not insulating 
but conductive, which was also mentioned elsewhere [229, 230, ^3l[. The effective 
permittivity of the mixture was expressed as in the form of Eq. (g) ll9|, with 



Ae 



£i[V£2 + (1 ~ v)£i] - g(l - v){£2 - £i) 
[r/e2 + (1 - ?7)ei] - m{£2 - £i) 

f]q{l - q){e2(Ji - £i(J2) 



[(1 - r])ei + 7?£2 -I- r]q{ei - £2)] [(1 - r])ai + 770-2 + Vli'^i ~ '^2W 
(1 - r])ei + ri£2 + r]q{ei - £2) 



(29) 



=£0 



(1 - f])ai + 770-2 + 779(0-1 - 0-2) 



where 77 is the inverse shape factor, rj — n^^. Note that for 7/ = 1 the equations above 
yield epxressions for Maxwell's laminated structure, Eq. (^ij). Similarly, Wagner's [ p6{ 
expression can be obtained by considering rj — 1/3. Steeman and Maurer have 
also introduced a third phase to the system where the third phase is an outer shell of 
the second phase (phase 3). 

The mean field approximations are valid for low concentration of inclusions 
{q <^ 1), since in such cases the inclusion particles are apart from each other and 
can not f eel the influence of the closest neighbors - dipole-dipole interactions [232, 
m, m, m, m, m, m, m, which accordingly are ignored. 



3.2.3. Molecular approaches In these approaches, the dielectric behavior of the 
mixture is assumed to be the summation of the dipole moments of each molecule 
in a vacuum matrix. Similarly to the mean field approach each molecule feels only the 
force of the applied field, not the neighboring molecules - no dipole-dipole interactions 
are included. Examples to the molecular approaches are Clausius-Mosotti [7| |24l|[ , 
Onsager-Bottcher ]242[ |37| , |^, |24l| , and the effective permittivities are expressed as 
in Eqs. ( p^ and (|29|). Kirkwood [243[ extended the Onsager theory of dielectric 
polarization by including the influence of nearest neighbors. 
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3.2.4- Regular arrangement of inclusions Rayleigh |244| approached to the binary 
mixture problem from a different point of view by assuming unit cells. In his 
approach, repeating unit cells composed of a matrix phase with sphere inclusions 
in the centers were considered. He assumed that both media were dielectrics without 
any conductivity, then, the permittivity of the composite was expressed analytically 



as a series: 



£l 



£2 



e2 + 2ei-g(e2-£i)-gi°/3^ 



6(52+45i/3) 



(30) 



Rayleigh's expression has also been derived by Lorentz-Lorentz 226 1 and Maxwell- 
Garnett [245|. The model includes the interaction between spherical inclusions, 
therefore, by contrast to Wagner ||l6| and Sillars |l^, it should be applicable to 
mixtures with higher concentration of inclusions. However, when the concentration of 
the inclusions is approaching one, q — > 1, the permittivity of the mixture, eis according 
to Eq. ( |3C| ) not equal to the permittivity of the inclusions, £2- Although the model 
can be used for irregular distributions of dil ute mixt ures, for highe r con centrations of 
inclusions this formula deviates from reality [246, 247 1. Bruggeman |248| used Eq. (^o|), 
and introduced the incremented volume fraction, dq, of the inclusion phase. He has 
then derived the following expression: 



(1-9) 



£1 



£1 + 2ee 



£2 



£2 + 2ee 







(31) 



This equation is known as the symmetric Bruggeman formula. The non-symmetric 
Bruggeman formula is p4^, |24g|]. 



£2 



£2 - £1 



(1-9) 



(32) 



Here again 77 can be regarded as the inverse shape factor as in Eq. ( p7| ) or space 
dimensionality |0, |^ . 

At higher filler concentrations the mutual interaction of particles is important. 
Interactions between two particles with different diameters and different dielectric 
properties which were placed in a background medium were presented by Emets [233| 
in two-dimensions. The mutual influence of neighboring particles, e.g., force, F, not 
only depends on the mentioned properties, but it is also affected by the direction of the 
applied electric field. The polarization of each particle in the applied field affec ts th e 
field distribution, and therefore, the polarization of the neighboring inclusions [239|. 

The Rayleigh model was improved by Emets |250| for systems with three 
components which had complex dielectric permittivites, ^ e^. In this case, the 
inclusions were assumed to be doubly periodic in an alternating order, and these were 
two different cylinders with radii, ri and r2. The complex effective dielectric constant 
of the system, e, was calculated from the averages of the displacement and electric 
field vectors in the region of interest. 



e =£1 



A 
B 



A=l- Ai2qi/2 - Ai3g2/2 + A^^Ai + Al^A2 + Ai2Ai3(Bi 



B=l + Ai2<zi/2 + Ai3(Z2/2 + Ai^Ai 



A23A2 -f Ai2Ai3(Bi 



B2) 
B2) 



(33) 
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Here 



£l — £p 



(-1 < Alp < 1) 



(34) 



where, (71.2 are the concentrations of the mclusion phases, (0 < qi^2 < !)■ The 
parameters and Bk (fc = 1, 2) are functions of the radii of the inclusions. 



Ak ^2rl\2rlY: 



171 — 1 

Tk + 2m 



00 OC 

EE 

n— 1 m— 1 

Tfe — 2m + 1 



rk — 2m 



{rk - 2m)2 - An'^ 



(rfc + 2m)2-4n2 (rfe - 2m + 1)2 - (2n - 1)2 

Tfc + 2m — 1 
{rk + 2m - 1)2 - (2n - 1)2 



(35) 



00 00 



Sfc =2r2_J2r3^ 

I T?J.= 1 

Tfc + 2m 



— rt - (2m - 

m— 1 ^ 



m - 1)4 + E E 



Tfc — 2m 



(rfe - 2m)2 - (n - 1)2 



7'fc — 2m + 1 



+ 



(rfe + 2m)2 - (?i - 1)2 (rfe - 2m + 1)2 - 4n2 

rfe + 2m — 1 
(rfe + 2m - 1)2 - 4n2 



(36) 



Ak and Bk values converge quickly to a constant level for n = m > 10. When 
the Rayleigh conditions are apphed, Eq. ( |3^ ) becomes equal to Eq. (pO|). Although 
these equations for regular structures with round shapes (spheres or cylinders) 
are valid for low concentrations, at high frequencies the sizes of the particles are 

|25£ 



playing an important role |g5l|, g5^, g5^, gO^, g5J, |25^, |25§1. Finally, other 
periodic structures have also been presented in the literature such as spherical 
247 , ^46|, rH , c ylindrical inclusions |25£, [250[ tili ngs (checkerboard) 
|26l( |262| p63| p64| and layered structures |2"65|, [266| [267 |. 



particles | l257|, |258| 
structures]|2l|,]260 



3.2.5. The spectral function approach A theoretical advance to the modeling of 
dielectric propertie s of binary composites ha s be e n d e veloped i ndependently by 
' I, |26|, |7^, |7l| and Milton |7|, |5§, |7|, |7|, |7|. The method is 



Bergman | 

called the spectral theory, in which a function of composite has been introduced as 
a compact way of representing data over a range of frequencies and it highlights the 
role of geometry in determining the effective properties |276, 277, 275]. In this model, 
the dielectric permittivity of the composite Eg is expressed as a function of dielectric 
permittivities of the constituents £1 and £2 and the geometry of the composite. 



3. 3. Bounds on electrical properties 

Difficulties in calculating the effective electrical properties of mixtures lead to int erest 
in obtaining bounds on these parameters. Wiener |27£| and Hashin-Shtrikman [ 280 | 
proposed bounds for two component systems. Wiener assumed layered structures 
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which had topological configurations parallel or perpendicular to the applied field 
direction. According to it 



£i£2[q£i + (1 - q)£2] < Ee < (1 - q)£i + q£2 



(37) 



These bounds are called absolute bounds meaning that the effective property, £e can 
not be present outside the region given in Eq. . If one of the phases forms inclusions 
(phase 2) in the other (ph ase 1 ), then, these two expressions can further be put together 
by introducing a factor n [^4l[ , which can be treated as the shape factor of the inclusion 
phase. 



£e - £l 



= q 



£2 - £l 



+ {n-l)ei £2 + (»^ - l)£l 



(38) 



Eq. (|^ reduces to Wiener bounds as expressed in E q. (^ for n = 1 and for n = cx) 



Moreover, Eq. ( pSh i s also the same as the Rayleigh |244| and Maxwell-Garnett [245| 
expressions, Eq. (|3C|). Landauer js) has mentioned that n is the space dimensionality, 
where n = 2 or n = 3. His approach may be valid, as in other mixture formulas, i.e., 
the Sillars derivation using the effective medium th eory for a composite with arbitrary 
shapes ||l^ or the Bruggeman assymetrical theory |24^ using the differential effective 
medium theory, the n-value represents the shape of the mono-dispersed inclusion 
particles. 

Hashin-Shtrikman bounds are narrower than Wiener bounds and they are derived 
for models of homogeneous and isotropic mixtures. In two-dimensions, the Hashin- 
Shtrikman bounds are expressed as 



£l + 



(£2-?i)-i + (l-g)(2?i) 



— < £ < £2 



1-q 



(£1 



£2 



■9(2£2)- 



(39) 



Other bounds have been proposed by Bergman |281 



27C|, |269[, Milton g73|, g72| 



and G olden and his collaborators [282, 283, 284, 285 1. It was found that (see Tuncer 
et al. |286|| ) the bounds were not valid for low frequency permittivity values when 
interfacial polarization was considered. 



3.4. Percolation theory 



When higher inc lusion concentrations or higher ratios of electrical properties of 
are considered, percolation phenomenon should be taken into 
288|, |28g|, |290|, |29l|, It is therefore more meaningful to speak of 



constituents |264 



account |287 



percolation in disordered systems 1 293 1 . At some specific concentration, depending on 
the model of the packing density, particles start to form conducting chains (percolating 
clusters) which infiuence the electric properties, the overall conductivity and the 
dielectric relaxation |287|. In regularly distributed systems composed of hard spheres 
(in 3-dimensions) or hard disk s (in 2-dimensions) percolation is only observed close 
to limiting concentrations |l^, 246, |2^, Other studies have shown that not only 



the percolation of one of the phases but the field, charge 
particle related issues [212, 192, 296, 297, BS, 298[ should also be considered 



239, 294, m, 295 [ and 



In the case of disordered systems, their physical behavior can usually be 
described by a power law with a percolation threshold and a critical exponent 
depending on the geometry [299, 300, 20£, |301[ and on the conduction of the 



system [287, 288, 291, 290, 302 1. One should also not exclude the local field changes 



Dielectric mixtures 



14 



in the composite due to interaction of inclusions. Sorbella [ |303[ concluded that the 
conductivity of the system around the localized scatterers (inclusions) did not obey 
Ohm's law. This shows that the micro-structure of a heterogeneous medium affects 
the conductivity of the system depending on a length scale, Ic which was defined 
as the distance between the scatterers or inclusions. Mendelson and Schwacher 
have also mentioned a similar length scale, Id, which is the variation of the dielectric 
permittivity in a continuous random dielectric. They have emphasized an other length 
scale, Lqi variation of charge density or the electric field. Therefore, average values 
for the physical quantities such as electric field, potential, dielectric permittivity etc. 
can be used when Lg ^ or Lq ^ Ic- 

High concentration of inclusi ons i n compos i te st r uctu res brings forth the problem 
of particle-particle interactions 



304 




The question of packing 
density for the particulated composite materials | 305| , 306 1 and of the packing structure 
itself becomes important. Particle size, shape, size distribution and cohesion force 
between them affect the maximum filler content |^55| and isotropy of the material. In 
the percolation region (concentration) and over, the theories described in this section 
become unreliable. However, some of the mixture formula might be used to predict 
the percolation threshold, in such cases the shape parameter n in the mixture formulas 
may be connected with the threshold . 

The electric conductivity and its mechanisms in granular materials containing 
metallic inclusions have been investigated intensively [294, 209, 206, 307, 308, 
Cases when the concentration of the metal inclusions were low enough 
206l 207, 30? I due to possible technological applications, 
easy to estimate metal concentration for the metal- 



208|. Cases when the 
were interesting [ [307| 
However, it has not 



insulator transition 



been 
09| 



290 



311 1, since the tunneling probabilities between 
the metal inclusions (or the more conductive phase) are involved | p94 |312{ . For low 
concenrations of inclusions, in which the system is under the percolation threshold, 
the conductivity can be activated by the tunneling process. When the percolation 
threshold is passed, the system clearly shows dc conductivity which is similar to 
transition between hopping and dc conductivity 



4. Numerical simulations 



Previous approaches 

Numerical simulations on the behavior of composites can be performed using different 
techniques. One of the ways is to solve Maxwell's equations in a computational 
domain |313, 314, 315| , ^16[ and to obtain either the electrical potential or electric field 
distributions in the domain. However, differences in the physical dimensions in the 
domain should not be too large, which is one of the drawbacks of computer simulations. 
Once the distribution of these quantities ar e kn own the effective material parameters 
can be calculated using different methods |317[ . The finite element method (fem) is 
one of the most intensively used approaches in the field simulations 313 314 , 318 1. In 
this study, the ef fecti ve medium p r operties of mi xtu res a re calculated using the finite 
element method ^l^, ^ |28^, [3191, ^ |32l|, |32^, ^ |32^, Si milar approaches can 

|32? 



also be found elsewhere ||32|, |30|. 




23L 324, 325 



233 



326 



232 



327 



23C, 231 



329|, |330|, |33|, g32| ^| 



204 



There is a multipole expansion approach |25C 
in which the considered inclusion phase is replaced by multipole source distributions. 
Monte Carlo techniques have been introduced in which many possible geometrical 
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arrangements have been considered |334, 335, 336 1. The Monte Carlo approaches are 
similar to molecular dynamic simulations. Since each medium can be characterized 
with its resitance, capacitance and inductance, modeling composite media with 
resistor-capacitor networks has studied g |33|, |3|, |39[ |34|, |o|, |34ll|342|. 
In addition, if there exists a periodicity in the considered geometries, i.e., arrangement 
of inclusions, fast-fourier-transformation techniques has been applied to obtain the 
response of the mixtures g, |l34[ |l3|, |l3^, |34|, |4|, |4|, |346[ p7|, Finally, novel 
approaches, like expanding the po tenti a l in ha.rmo n ic hmction s |34£ , 35C| , pi , 351, 352 1, 



using Gaussian random fields 



353, 354 



362, 363 



355 



15C, 364, 365, 142, 



356, 357, 345 



and 



m any 



other 



366| , |367| , |368i have 



methods [g5§ g5§ gSg g60l g61 
been introduced. 

At high frequencies or in structures where the smallest part of inhomogeneities - 
which can be regared as the size of inclusions - are larger than the wavelength of the 
applied field, the effective permittivity calculations are not trivial and the results can 
be misleading [251, ^69, 



Numerous simulations confirmed the validit y of the mixtu re formul as for reg ular 
structures when the mixtures were dilute |5 



21 322, 323 



3 324 



231, 327, 363 1 . For dilute random mixture the results were also similar |371 



32, 
171 

emphasizing the validity of the effective medium approximations at low inclusion 
concentrations. When higher concentration of inclusions were studied the validity 
of the analytical solutions and their comparisons to numerical results were usually 
questioned and the results were compared with those of the bounds on the 
electromagnetic properties of composites |27£, 



281, 274, 273, UTA 27C, 271 



28C, 282, B72 283, 284, 285, 373 



691 374, 375|. In addition, at higher concentrations, 
approaches like fundamental geometrical parameters have been introduced [272, 275, 
|345|, 281, 27C, 271 [. Although the geometrical parameters of a structure were 



constant, one should not disregard the effect of the ratio of electrical properties of 
the constituens [p4l|, 342, 376 1 as mentioned previously. Finally, in the numerical 



calculation, little attention has been given to the interfacial problem |377, 378 



4.2. Application of the finite element method 

Analytical calculations of electromagnetic problems are limited to geometrical 
constraints. For some simple geometries with a small number of materials (regions) 
and symmetries, analytical solutions can be found [229, 233, 25C, |379[ ^T\ . These 
analytical solutions are obtained using methods of images |380, 381 1, orthogonal 
functions (Green functi ons) [ |382{ and complex variable techniques (conformal 



mapping) [(3831 , |381 



229[. For more complex geometries and for non-homogeneous 
regions composed of several materials, numerical solu tions of partial differential 
equations and of integral equations have been developed [^313) . 

Numerical solutions of electrostatic problems within a non-conducting medium 
are based on solving Poisson's equation 



V • (eeoV0) 



(40) 



where </> and p denote the electrical potential and the charge density in the considered 
region, respectively. If the medium is conductive where no free charges and sources of 
charges are allowed, then, the solution is given by 



V • (crV0) = 



(41) 
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n = 3 



n = 4 




• • • 




Figure 4. Finite-size scaling of the inclusion shapes (polygons). 



When the medium is a mixture of these two cases (lossy dielectric), it consists of 
dielectric and conductive components. The solution, then, becomes time dependent 
and is given by a complex electric potential in the region with the coupling of Eqs. 
( ^o|) and (|4l|), which is also known as the continuity equation for the current density 



V-(aV0) + ^{V-[(e£oV0)]}-O 



(42) 



Equivalently in Fourier-space with frequency dependent properties, i.e., the complex 
dielectric permittivity, e{ijj) = e + x(a;) -I- (j/{ieQUj), where no free charges are allowed 
in the region, due to conductivity of the medium (lossy dielectric) 



V • {[ieoe{uj)uj] V0} = 



(43) 



Field calculation softwares, such as ACE and Femlab 384 , 385 , can be used to solve 
Eq. (^^. Once the potential distribution is known, the complex permittivity of a 
heterogeneous medium can be calculated in several ways, e.g.^ 

(1) by using the total current density, j, and the phase difference, [386, 3^, |5, 286 |, 

(2) Gauss' law and losses [|1|2| |||86|, 



(3) energy balance |286 



(4) and by using the average values of dielectric displacement (D) and electric field 
(E) g,|3,|28§. 

To check the accuracy of the finite element results, two-dimensional inclusions in 
the form of re(7uZor polygons with n sides were considered [322|, as displayed in Fig. ^. 
The phase parameters, e and a, were frequency and voltage independent. The values 
of e and a were chosen such that the interfacial polarization observed had a relaxation 
time, T, around 1 s. This was achieved when the matrix phase had ei = ei = 2 and 
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Figure 5. Dependence of percentage error in liigli frequency dielectric 
permittivity, |£m — A|/A on number of regular polygon sides, n. The symbols 
open (O) S'lid filled (•) indicate tlie solutions obtained using the current density 
and phase shift between the applied voltage and current and using Gauss' law 

and the total losses in the medium, respectively. The solid lines ( ) represents 

the fitted curves, ric is the critical side number for regular polygons. 



(Ti = 1 pS/m, and the inclusion phase had 62 = £2 — 10 and (T2 = 100 pS/m. The 
polygons were generated using a circle with radius, r, and a constraint on the area of 
the polygons, q. Then, the radius, r, as a function of n and q is expressed as. 



The denominator inside the square root approaches tt as n — > cxd. The size of radius 
was also used for the meshing procedure of the computation domain where r„/15 
was the size of the minimum triangle. Furthermore, this approach led to a finite-size 
scaling |5^, |387| in which as n — > 00, the inclusion phase was a perfect disk. In Fig. ^, 
an example of the finite-size behavior is shown. In the figure, difference between the 
dielectric constant at high frequencies Em = s{uj@lMHz) and the value obtained from 
analytical solution in Eq. (|3|) which is presented as A, (A = e^g. Q(w@lMHz)), is 
plotted against the number of sides, n. A critical number of sides, ric, is defined, 
such that over this value, n > ric, the effective properties of medium with regular 
polygons as inclusions are approximately similar to those of a medium with disk 
shape inclusions. The analysis showed that is approximately 15 regardless the 
concentration levels considered, q < 0.3, and the error in the calculations is < 0.01 % 
for Uc > 15, as displayed in Fig. ^. In fact, even a decagon {n — 10) can imitate a 
disk in which the error in the calculated electrical quantities is less than < 0.1 %. 




(44) 
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Figure 6. Comparison of tiie calculated normalized dielectric responses of a 
binary composites for q = [0.1(*), 0.4(»), 0.7(H)] and the analytical solutions for 
q = 0.4. Properties of the media are ei = 2, £2 = 10, a"i = lO"-'^'^ Q~^m~^ and 

(72 = 10"^" il~^m~^. The solid line ( ) represents a Debye-type dielectric 

response. 



4-3. Concentration of inclusions and dielectric relaxation 

By scanning the dielectric properties of mixtures in a frequency window, on gets a 
better understanding of the relaxation and the importance of the local electric field or 
in other words the interaction between inclusions. Therefore, analysis and comparison 
of mixtures with intermediate (g = 0.4) and high (g = 0.7) concentrations of inclusions 
were discussed in Ref. |]55|] . The inclusions were hard disks which did not overlap; this 
resulted a packing density (or limiting concentration) for the model, which was 7r/4 
for the two dimensional square lattice. The obtained frequency dependent complex 
dielectric permittivities were compared to analytical mixture formulas and were also 



curve-fitted [388| with the addition of two contributions; the Cole-Cole empirical 

(45) 



expression [[76[ and the dc conductivity, 

X(0) , (a/eo) 



1 -I- {iWt)"' lUJ 

The main reason for the curve-fitting procedure was to check the behavior of dielectric 
relaxation and to find a plausible explanation for deviations from exponential dielectric 
relaxation. The fitting procedure is presented in Ref. |Q . 

4. 3.1. Intermediate concentrations of inclusions When the concentration of the 
inclusions phase was 0.4, the differences between different mixture formulas and the 
FEM solution were significant, as displayed in Fig. ^ and the curve fitting parameters 
are presented in Table |l|. In the figure, the dc contributions have been substructed from 
the data. The analytical solution of Emets |25C| ] (Eqs. ( p3| ) and (p^)) was separated 
from the other formulas. There were no differences between the responses obtained 



from the Wiener |279 and Steeman |l9] expressions. Except for the Bruggeman [248| 



formula, the other dielectric responses had similar high frequency values, e, however. 



the Emets [250| solution yielded the highest value of £^@lMHz) compared to the 
others. The fem solution was close to those Wiener ^79[ and Steeman (l9| but did 
not have exactly the same shape, compare the fitting parameters in Tablelu 



Dielectric mixtures 



19 



Table 1. Cole-Cole fitting parameters of different dielectric mixture formulas and 
the finite element calculations. The electrical parameters of the media are ei = 2, 
£2 = 10, CTi = 1 pSm-l and (T2 = 100 pSm"!; q = 0.4. 



Model 



x(0) 



a/eo 



Steeman 
Bruggeman [248| 
Emets 125^ 



Wiener 12791] 

FEM 



3.4545 
3.5785 
3.4881 
3.4545 
3.4547 



1.0447 
1.6213 
1.1589 
1.0447 
1.0584 



1.2689 
1.6331 
1.3305 
1.2687 
1.2658 



1.0003 
0.9988 
0.9993 
1.0003 
0.9897 



0.2586 
0.3031 
0.2677 
0.2586 
0.2594 



The relaxation times, r, from the curve fitting procedure were close to each other 
except for the Bruggeman [p48[ response. The conductivity terms, a/eo obtained from 



the curve-fitting procedure also showed that the Bruggeman ]24^ and Emets [250| 
formulas yielded higher values than the others. 

4-3.2. High concentration of inclusions When the concentration of the inclusions was 
0.7 which is close t o th e limiting concentration of the square lattice, q — tt/A — 0.7854, 



only the Wiener |27£] and Steeman |l9[ responses were similar to each other, as 
presented in Table 2|. All the responses were different, and the high frequency dielectric 
permittivities e were between 5.5 and 6. The low frequency values of permittivity, 
£^(a;@l[j.Hz), on the other hand, had large discrepancies. It was extraordinary that 



the Emets |25C| formula resulted in the highest value of the real part of the dielectric 



permittivity at low frequencies, e -I- x(0), while in the other considered cases the 



response estimated by Bruggeman [248] had this character. The fem solution was 
non-symmetric, and it could not be modeled with the Cole-Cole empirical formula as 
presented in Fig. |[ The fitting parameters of the fem have shown that the obtained 
reponse did not have any Cole-Cole characteristics. 

Thus, for a simple computational domain of a binary structure, the dielectric 
relaxation was Debye-like when the concentration of the inclusions was low, and the 



mixture formulas of Steeman Jh9l, Wiener |27£] and Emets ]E50l agreed with the fem 



solutions |5^. Non-Debye and non-Cole-Cole dielectric relaxations were observed for 
the FEM solutions at higher inclusion concentrations. Since a simple binary system 
composed of a square lattice exhibited non-expected dielectric responses, introducing 
geometries and structures which remind more realistic situations would probably yield 
non-Debye like responses as well. When the dimensions and shapes of the inclusions 
were different, the time constant of the relaxation would not be single but distributed. 



4.4- Structural differences and dielectric relaxation 

4.4- 1- Ordered structures with different lattices The considered regular lattices were 
square and triangular as presented in Fig. ^ In these structures, the smallest repeating 
units used in the calculations are shown as shaded areas in the figure. Two different 
cases, in which the inclusions were more (cti < a2 - match-composite) and less 
((Ti > (72 ~ reciprocal-composite) conductive than the matrix phase were taken into 
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Table 2. Cole-Cole fitting parameters of different dielectric mixture formulas and 
the finite element calculations. The electrical parameters of the media are ei = 2, 
e2 = 10, CTi = 10-12 n-im-l and (T2 = 10'^° Q-^m-^; q = 0.7. 



Model 



X(0) 



cr/eo 



Steeman 
Bruggeman [248| 
Emets 125^ 



Wiener 12791] 

FEM 



5.4998 
5.8779 
6.0313 
5.4998 
5.7991 



4.7154 
10.056 
10.624 
4.7153 
6.9035 



1.7875 
3.4816 
3.1293 
1.7875 
2.1502 



1.0003 
0.9651 
0.9975 
1.0003 
0.9609 



0.6068 
1.0528 
1.0697 
0.6068 
0.7685 



• •• 




Figure 7. Square (left) and triangular (right) lattices with square and hexagonal 
unit cells. Dark circles are the inclusions which are distributed in the host medium 
(white background). The shaded areas are the primitive (unit) cells which are used 
in the calculations. 



consideration. The ohmic conductivities of the media were either 1 pS/m or 100 pS/m. 

The analysis of the dielectric responses by means of the Monte Carlo technique 
are presented in Fig. ||a. The product of dielectric strength and the distribution of 
relaxation times, Ae x JT, of the ai < case, (0 — □), for two different regular 
structures were similar and nearly-symmetrical. The slopes of the distributions could 
be fitted by (logr)'' expression, where k 6 for logr < 0, and k « —5 for logr > 0.3. 
It is possible to convolute this distribution into two separate ones characterized by 
different dominant time constants. When txi > cr2, (• — ■), the distributions were 
broader. Moreover, the slope of the distributions at shorter r-values were identical 
but, the arrangements of the inclusions had altered the distribution of relaxation times 
at longer times, logr > —0.1. The reconstructed dielectric responses are presented in 
Fig. ||b and they were more successful than the Cole-Cole curve-fittings. 



^.^.2. Disordered structures Earlier numer ical simulations^, pi |^ ||, |32|, |46[ ||, 



134 1 of dielectric mixtures by means of the fem |22, 20l p5, 325, 23 have not considered 



frequency dependent properties of disordered systems. However, some experimental 
and theoretical investigations have been reported |287, ^, 339, 337, 340 1. The authors 
of this report, as in regular lattice structures investigated behavior of random systems 
with concentrations of the phases constant, gi,2 ~ 0.5. A concept by of puzzle-like 
structures was introduced, as presented in Fig. |^. The networks were composed of 
16 X 16 cells, and the matrix and inclusion elements were placed randomly. Altogether 
1024 structures were analyzed. A range of responses were obtained depending on the 
ratio of cunductivity values of the structure constituents as well as their distribution. 
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Figure 8. (Up) Product of dielectric strengtli, Af , a nd distribution of relaxation 
times, jr, found using the Monte Carlo technique p2\ . The open (o — □) and filled 
(• — I) symbols represent the responses when ai < cr2 (match-composite) and 
fi > (72 (reciprocal-composite), respectively. The boxes (□— ■) and circles (o— □) 
represent the square and triangular lattice, respectively. (Down) Reconstructed 
dielectric responses, the chain lines ( — ■ — ), using the distributions in upper 
subfigure. 



Examples of the frequency dependent dielectric permittivities of the structures 
are displayed normalized in Fig. |l^ after subtracting of the ohmic conductivity 
contributions. The effect of randomness was significant compared to the regular 
structures, the influence of the local fleld enhancement and internal structure yielded 
stronger polarizations. The curve fitting parameters given in Table ^ have indicated 
that there were slight differences between the two empirical formulas for the txi < (J2 
case (Fig. ^0|a), however for the cti > (J2 case (Fig. |10|d) the Davidson-Cole expression 
produced a more successful fit (lower residual) than the Cole-Cole one. 

All the obtained responses started and ended at the e'-axis with an angle of 7r/2, 
which indicated that they were neither Cole-Cole nor Davidson-Cole; a Cole-Cole 
response starts with aw/ 2 and, similarly, a Davidson-Cole response starts with (37r/2 
in the Cole-Cole representation. This has shown that there exist two relaxation time 
constants, r, which set boundaries for the distribution of relaxation times | |5l[ |. The 
reconstructed response using the distribution of relaxation times is displayed by lines 
in Fig. I^. 
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Figure 9. Structures in which the minimum, e'(ui@2007r Hz), (a&c) and the 
maximum, e'{w@2iT mHz), (b&d) dielectric permittivities are obtained when 
o"i < CT2 (a&b) and ai > cr2 (c&d). 



Table 3. Curve fitting parameters of Eq. ( psj ) using Cole-Cole and Davidson- 
Cole empirical expressions for the dielectric relaxation of random structures with 
square networks. 



Structure case e Ae r[s] a a/eo[S/F] 

Fig. |a 0-1 < (72 4.145 4.819 3.916 0.921 0.609 

Fig. 9b 0-1 < (72 4.654 9.260 6.303 0.927 1.185 

Fig. 9c (71 > (72 4.424 33.89 9.329 0.751 1.623 

Fig. |d (71 >(72 4.803 65.29 23.49 0.818 1.159 

^ 

Fig. |a (7i<(72 4.134 4.796 5.045 0.818 0.609 

Fig. 9b (71 < (72 4.630 9.222 8.068 0.822 1.185 

Fig. 9c (71 > (72 4.085 33.44 23.88 0.523 1.623 

Fig. |d (71 >(72 4.213 64.66 46.62 0.600 1.159 



The analysis has showed that the randomness in the medium spread the dielectric 
relaxation times over a broder region compared to the regular structures. If the finite- 
size effects had been taken into account, assuming that the calculations were done 
on larger crosswords (not 16 x 16), or if our results could have been extrapolated, 
the dielectric relaxation in random media would yield a behavior reminding the low 
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frequency dispersion ^89| , p8[ . 

4- 4- 3. Micro structural differences in random lattices When random structures with 
different lattices, as shown in Fig. were considered, square and honeycomb tiles 
the smallest parts. The frequency dependent complex dielectric susceptibilities of 
selected structures with different lattices having the same distribution of phases are 
presented in Fig. [T^. The shape of the Cole-Cole plots of the obtained responses 
after the substructing of the ohmic losses were non-exponential, similar to those 
presented previously in Fig. ^ as presented in Fig. |l2|. In the figure, the reponses 
were for the reciprocal-composites (cti > 02), and in the investigation two different 
conductivity ratios of the constituents were taken into account. Only the structure 
with square tiles had symmetrical response for ei/£2 = 1/5 and oxjai — 100, 
the other reponses were asymmetrical. When the ratio of the conductivities were 
increased, ai/(J2 = 1000, the dielectric relaxation of structure with square tiles 
became distorted. It was interesting that the structure with honeycomb tiles had 
similar dielectric responses for two considered cases with different conductivity ratios. 
When the conductivity ratio was increases, a condition for the 'normal' low frequency 



dispersion ||3^, 337, 389 was created since both responses had flat regions. The ohmic 
conductivities of the structures were different and the structure with honeycomb tiles 
had higher conductivity values. 

5. Discussion and future challenges 

In dilute systems, the dielectric relaxation due to interfacial polarization yields 
classical Debye relaxation, which was also verified by numerical simulations. However 
for dense mixtures non-Debye (non-exponential) relaxation behavior can be observed 
depending on the internal structure of the system and the electrical properties of the 
constrituents. The differences between regular and random structure were remarkable, 
such that, the random structures had broader non-Debye dielectric relaxations when 
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Figure 10. Normalized dielectric responses of random structures with square 
network. The ininimum and maximum dielectric permittivities obtained when 
fl < (T2 (Fig. ^ & b) and ui > 02 (Fig. gc & d). The opeii (o — □) and 
filled (• — H) symbols represent the responses of structures in Fig. (3 in which the 
lowest, £(u'@2007r Hz), and the highest, £{w&2tt mHz), dielectric permittivities 
are obtained, respectively. 
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Figure 12. Normalized dielectric responses obtained for two chosen random 
distributions with square (□ — ■) and honeycomb (A — A) tiles. The open (□— A) 
and filled (I — A) symbols represent the real, x' a^nd imaginary, x" parts of the 
calculated dielectric susceptibility. 



the complex dielectric permittivity ratio of the phases were small. Further modification 
of the constituents' properties indicated that, it was even possible to obtain a dielectric 
response remanding the low frequency dispersion behavior. Therefore, one should be 
careful in interpretating dielectric data since the interfacial polarization could result 
in various dielectric relaxations depending on the internal structure of the mixture 
and on the properties of the constituents. 

The most important question arising which concerns any kind of numerical 
simulations is the reliability of obtained results. Here the reliability does not mean 
the correctness of a physical model representing a real physical phenomenon, but 
its mathematical and computer implementations. In the case of modeling dielectric 
properties of composites, this is one of the basic problems because usually there are 
no exact analytical solutions available that would allow for checking the accuracy of 
the computed data. Special attention should be paid, for example, when a method 
for calculating the effective complex permitivity is chosen. Four different methods 



mentioned in § 4.2 were compared, but only two of them (one based on the total 



current and phase shift between the current and applied voltage; and one where 
average field and dielectric displacement were calculated) yielded stable and reliable 
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solutions J317| , |322{ . Physical (or mathematical) reasons of discrepancies of the data 
obtained with different methods are not completely understood yet. 

It is obvious, that correct numerical representation of the initial problem must 
be used in simulations and appropriate meshing of the computational domain must 
be performed. Meshing is extremely important for the considered problem in order to 
resolve anisotropy in the composite structure. In this connection, special care should 
be taken about meshing ability of the used software for solving the field problem 
in Eq. (p3|), because it could be very difficult and even impossible to generate a 
computational mesh if great differences in the dimensions of the matrix and inclusions 
exist. 

The desired outcome from simulations is an agreement between computed results 
and characteristics of a real composite material. That would prove the applicability 
of the model and allows performing simulations for other structures with unknown 
effective parameters. This can be achieved if a real three-dimensional structure of the 
material is built, and the dependences of electrical properties of the constituents on 
state variables (frequency, temperature, pressure, electric field, etc.) are introduced in 
simulations. Furthermore, a possibility of modeling different properties of interfacial 
layers between constituents of the composite should be introduced. These problems 
have not been solved yet, and is the subject of on-going research. 
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